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States of self-stress, tensions and compressions 
of structural elements that result in zero net 
forces, play an important role in determining the 
load-bearing ability of structures ranging from 
bridges to metamaterials with tunable mechanical 
properties. We exploit a class of recently intro¬ 
duced states of self-stress analogous to topolog¬ 
ical quantum states to sculpt localized buckling 
regions in the interior of periodic cellular meta¬ 
materials. Although the topological states of self 
stress arise in the linear response of an idealized 
mechanical frame of harmonic springs connected 
by freely-hinged joints, they leave a distinct sig¬ 
nature in the nonlinear buckling behaviour of a 
cellular material built out of elastic beams with 
rigid joints. The salient feature of these localized 
buckling regions is that they are indistinguish¬ 
able from their surroundings as far as material 
parameters or connectivity of their constituent 
elements are concerned. Furthermore, they are 
robust against a wide range of structural pertur¬ 
bations. We demonstrate the effectiveness of this 
topological design through analytical and numer¬ 
ical calculations as well as buckling experiments 
performed on two- and three-dimensional meta¬ 
materials built out of stacked kagome lattices. 

Mechanical metamaterials are artificial structures 
whose unusual properties originate in the geometry of 
their constituents, rather than the specific material they 
are made of. Such structures can be designed to achieve 
a specific linear elastic response, like auxetic (negative 
Poisson ratio) [T] or pentamode (zero shear modulus) [2] 
elasticity. However, it is often their nonlinear behaviour 
that is exploited to engineer highly responsive materi¬ 
als, whose properties change drastically under applied 
stress or confinement PHZ]. Coordinated buckling of the 
building blocks of a metamaterial is a classic example of 
nonlinear behaviour that can be used to drive the aux¬ 
etic response 015 ], modify the phononic properties [5] or 
generate 3D micro/nanomaterials from 2D templates [U] . 

Buckling-like shape transitions in porous and cellular 
metamaterials involve large deformations from the initial 
shape, typically studied through finite element simula¬ 
tions. However, many aspects of the buckling behaviour 
can be successfully captured in an approximate descrip¬ 
tion of the structure that is easier to analyze [iiniiiQiin]. 
Here, we connect the mechanics of a cellular metamate¬ 
rial, a foam-like structure made out of slender flexible el¬ 
ements da US], to that of a frame —a simpler, idealized 
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assembly of rigid beams connected by free hinges—with 
the same beam geometry. We exploit the linear response 
of a recently introduced class of periodic frames da, in- 
spired by topologically protected quantum materials, to 
induce a robust nonlinear buckling response in selected 
regions of two- and three-dimensional cellular metamate¬ 
rials. 

Frames, also known as trusses, are ubiquitous minimal 
models of mechanical structures in civil engineering and 
materials science. Their static response to an external 
load is obtained by balancing the forces exerted on the 
freely-hinged nodes against internal stresses (tensions or 
compressions) of the beams. However, a unique set of 
equilibrium stresses may not always be found for any 
load [15]. First, the structure may have loads that can¬ 
not be carried because they excite hinge motions called 
zero modes that leave all beams unstressed. Second, the 
structure may support states of self-stress - combina¬ 
tions of tensions and compressions on the beams that 
result in zero net forces on each hinge. An arbitrary lin¬ 
ear combination of states of self-stress can be added to 
an internal stress configuration without disrupting static 
equilibrium, implying that degenerate stress solutions ex¬ 
ist for any load that can be carried. The respective counts 
Am and N^s of zero modes and states of self-stress in a 
d-dimensional frame with nodes and nb beams are 
related to each other by the generalized Maxwell rela¬ 
tion [mini: 

dn^ 77/b — Am Agg, (1) 

As Eq.[2 shows, states of self-stress count the excess con¬ 
straints imposed by the beams on the dn^ degrees of free¬ 
dom provided by the hinge positions. 

States of self-stress play a special role in the mechan¬ 
ical response of repetitive frames analyzed under peri¬ 
odic boundary conditions [18]. Macroscopic stresses in 
such systems correspond to boundary loads, which can 
only be balanced by states of self-stress involving beams 
that cross the boundaries. If these states span the en¬ 
tire system, boundary loads are borne by tensions and 
compressions of beams throughout the structure. Con¬ 
versely, by localizing states of self-stress to a small por¬ 
tion of a frame, load-bearing ability is conferred only to 
that region. Our approach consists of piling up states 
of self-stress in a specific region of a repetitive frame so 
that the beams participating in these states of self-stress 
are singled out to be compressed under a uniform load at 
the boundary. In a cellular material with the same beam 
geometry, these beams buckle when the compression ex¬ 
ceeds their buckling threshold. 

Although our strategy is of general applicability, it 
is particularly suited to isostatic lattices nnmui, with 
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dn^ = rib and no zero modes other than the d rigid body 
motions under periodic boundary conditions. According 
to Eq. these lattices only have d states of self-stress, 
insufficient to bear the d{d + l)/2 possible independent 
macroscopic stresses dam]- Additional states of self¬ 
stress can significantly modify the load-bearing ability 
of these structures at the cusp of elastic stability. Triv¬ 
ial states of self-stress can be created by locally adding 
extra beams between hinges. However, some isostatic pe¬ 
riodic lattices can harbour topological states of self-stress 
which owe their existence to the global structure of the 
frame. These frames, introduced by Kane and Luben- 
sky [H], are characterized by a polarization Rx whose 
origin can be traced to topological invariants calculated 
from the geometry of the unit cell. Domain walls between 
different topological polarizations as in Fig. la dill], 
and lattice defects [23], can harbour localized states of 
self-stress which can be used to drive localized buckling. 
Unlike their trivial counterparts, the existence of these 
topological states of self-stress cannot be discerned from 
a local count of the degrees of freedom or constraints in 
the region. An attractive feature for potential applica¬ 
tions is their topological protection from perturbations 
of the lattice or changes in material parameters that do 
not close the acoustic gap of the structure dlUHlZ]- 
In the remainder of this article, we use robust states of 
self-stress to design buckling regions in topological meta¬ 
materials composed of flexible beams rigidly connected 
to each other at junctions. As shown in Fig. 1, the states 
of self-stress are localized to a quasi-2D domain wall ob¬ 
tained by stacking multiple layers of a pattern based on a 
deformed kagome lattice HU- The domain wall separates 
regions with different orientations of the same repeating 
unit (and hence of the topological polarization Rx of the 
underlying frame). When the material is compressed uni- 
axially, beams participating in the topological states of 
self-stress will buckle out of their layers in the portion of 
the lattice highlighted in red, see Supplementary Movie 1. 
This region primed for buckling is indistinguishable from 
the remainder of the structure in terms of the density of 
beams per site and the characteristic beam slenderness. 
We verify using numerical calculations that the states of 
self-stress in the idealized frame influence the beam com¬ 
pressions of the cellular structure in response to in-plane 
compressions along the edges. The phenomenon survives 
even when the patterns on either side of the domain wall 
are nearly identical, reflecting the robustness of the topo¬ 
logical design. 


I. LINEAR RESPONSE OF THE FRAME 

The distinguishing feature of the isostatic periodic 
frame used in our design is the existence of a topological 
characterization of the underlying phonon band struc¬ 
ture m- If the only zero modes available to the structure 
are the d rigid-body translations, its phonon spectrum 
has an acoustic gap; i.e. all phonon modes have non- 
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FIG. I. Topological buckling zone in a cellular meta- 
material. a, Isostatic frame containing two domain walls 
that separate regions built out of opposite orientations of the 
same repeating unit (boxed; the zoom shows the three unique 
hinges in the unit cell as discs). This unit cell carries a polar¬ 
ization Rt = ai (solid arrow), and the periodic frame displays 
topological edge modes [T3]. The left domain wall harbours 
topological states of self-stress, one of which is visualized by 
thickened beams identifying equilibrium-maintaining tensions 
(red) and compressions (blue) with magnitudes proportional 
to the thickness. The right domain wall harbours zero modes, 
one of which is visualized by green arrows showing relative dis¬ 
placements that do not stretch or compress the beams. Modes 
were calculated using periodic boundary conditions, b, A 3D 
cellular metamaterial is obtained by stacking four copies of 
the beam geometry in a, and connecting equivalent points 
with vertical beams to obtain a structure with each interior 
points connecting 6 beams. The beams are rigidly connected 
to each other at the nodes and have a finite thickness. Points 
are perturbed by random amounts in the transverse direction 
and a small offset is applied to each layer to break up straight 
lines of beams. A 3D-printed realization of the design made 
of flexible plastic (see Materials and Methods) is shown in c. 
The sample has a unit cell size of 25 mm and beams with cir¬ 
cular cross-section of 2 mm diameter. The stacking creates a 
pile-up of states of self-stress in a quasi-2D region, highlighed 
by dotted lines. 
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zero frequencies except for the translational modes at 
zero wavevector. A gapped isostatic spectrum is charac¬ 
terized by d topological indices one for each primitive 
lattice vector {i G {1,2} for the 2D lattice used in our 
design). Although the phonon spectrum can be changed 
by smoothly deforming the unit cell geometry, the in¬ 
dices themselves are integer-valued and thus invariant to 
smooth perturbations. Their value can only be changed 
by a deformation which closes the acoustic gap. Analo¬ 
gously to topological electronic systems such as quantum 
Hall layers and topological insulators which harbour pro¬ 
tected edge states [28], the existence of these invariants 
guarantees the presence of zero modes or states of self¬ 
stress localized to lattice edges or domain walls where the 
rii change value. 

The count of topological mechanical states at an edge 
or domain wall is obtained via an electrostatic analogy. 
The lattice vector Rt = (Fig. la) can be in¬ 

terpreted as a polarization of net degrees of freedom in 
the unit cell. Just as Gauss’s law yields the net charge 
enclosed in a region from the flux of the electric polar¬ 
ization through its boundary, the net number of states of 
self-stress (minus the zero modes) in an arbitrary portion 
of an isostatic lattice is given by the flux of the topolog¬ 
ical polarization through its boundary M- In Fig. la, 
the left domain wall has a net outflux of polarization 
and harbours topological states of self-stress (the right 
domain wall, with a polarization influx, harbours zero 
modes). Although only one of each mechanical state is 
shown, the number of states of self-stress and zero modes 
is proportional to the length of the domain wall. Similar 
results can be obtained by using other isostatic lattices 
with an acoustic gap and a topological polarization, such 
as the deformed square lattice in Ref. |23| . 

The linear response of a frame can be calculated from 
its equilibrium matrix A, a linear operator which relates 
the beam tensions t (a vector with as many elements 
as the number of beams nb) to the resultant forces on 
the nodes p (a vector with one element for each of the 
2nn degrees of freedom of the rin nodes) via At = p. 
States of self-stress are vectors which satisfy At^ = 0 , 
i.e. they are beam stresses which do not result in net 
forces on any nodes (the index q identifies independent 
normalized states of self-stress that span the nullspace 
of A). The same null vectors are also an orthogonal 
set of incompatible strains of the structure, i.e. beam 
extensions that cannot be realized through any set of 
point displacements m- 

Since we are interested in triggering buckling through 
uniform loads which do not pick out any specific region 
of the lattice, we focus on the response to affine strains, 
where affine beam extensions ea are imposed geometri¬ 
cally by some uniform strain Cij across the sample via 

c.a=rTeijr^. ( 2 ) 

Here, a indexes the beams and is the end-to-end vec¬ 
tor of beam a. To attain equilibrium, the beams take 
on additional non-afhne extensions ena- Under periodic 


boundary conditions, affine strains are balanced by loads 
across the system boundary rather than loads on spe¬ 
cific nodes, which implies that the resultant beam ten¬ 
sions ta = k{eg, + Gna) must be constructed solely out of 
states of self-stress [HHS] (we assume for simplicity that 
all beams have identical spring constant k). Therefore, 
ta = = Tx whcrc T = [t^, P, ..., and the 

Xq are the weights of the N^s states of self-stress. These 
weights are determined by requiring that the non-afhne 
strains have zero overlap with the set of incompatible 
strains [ig (the affine strains are automatically compat¬ 
ible with the affine node displacements): 


f Te„a = f T Qtx - 

= 0, 

(3) 

which gives the solution 



X = /cT^Ga 


(4) 

^ta = fcTf'rea = fc^(t« 

• ea)t«. 

(5) 


q 


Therefore, the beam tensions under an affine deforma¬ 
tion are obtained by projecting the affine strains onto 
the space of states of self-stress. 

Eq. shows that the loading of beams under affine 
strains is completely determined by the states of self¬ 
stress. In a frame consisting of a single repeating unit 
cell, loads are borne uniformly across the structure. How¬ 
ever, if the structure also has additional states of self¬ 
stress with nonzero entries in conhned to a small re¬ 
gion of the frame, it would locally enhance tensions and 
compressions in response to affine strains, provided the 
have a nonzero overlap with the affine bond exten¬ 
sions imposed by the strain. Fig. 2 shows the approxi¬ 
mate states of self-stres^ for the domain wall geometry 
of Fig. la which have a nonzero overlap with affine exten¬ 
sions Gx and Gy due to uniaxial strains Sij = SixSjx and 
^ij = ^iy^jy respectively. The system-spanning states 
of self-stress and U in Figs. 2a-b respectively do not 
single out any particular region. Although they have a 
nonzero overlap with both Gx and Gy, they do not pro¬ 
vide significant stiffness to a uniaxial loading because 
a combined affine strain Sij = SixSjx + I36iy6jy exists 
with (3 ~ —1.2 such that the overlaps of the correspond¬ 
ing affine extensions with and U are small (less than 
10“^). Upon compression along one direction, say the 
frame can expand in the perpendicular direction to keep 
the tensions and compressions low in the majority of the 
sample. In contrast, the localized states of self-stress 
shown in Fig. 2c-d have a significant overlap with Gy 
but not Gx- Therefore, according to Eq. a uniform 
compression applied to the lattice along the vertical di¬ 
rection, with the horizontal direction free to respond by 


^ These approximate states of self-stress become exact when the 
separation between the two domain walls becomes large. 
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FIG. 2. Extended and localized states of self-stress of 
the frame, a—d, States of self-stress in the infinite periodic 
frame obtained by tiling the design of Fig. la in both direc¬ 
tions. The states a, b are largely uniform over the structure, 
whereas c, d are localized to the left domain wall. The over¬ 
laps of each state of self-stress with the affine strains Ox and 
Gy (see text) are also shown. Only states of self-stress with 
significant overlaps are shown; all other states of self-stress in 
the structure have • G{x,y} < 10~^. 


expanding, will significantly stretch or compress only the 
beams participating in the localized states of self-stress. 

Whereas the topological polarization guarantees the 
presence of states of self-stress localized to the left do¬ 
main wall, their overlap with one of the three indepen¬ 
dent affine strains is determined by the specific geometry 
of the hinges and bars. For the frame in Fig. 2, the 
states of self-stress visualized in Fig. 2c-d are crucial in 
triggering buckling response, which may be predicted to 
occur under compression along the y direction (or ex¬ 
tension along the x direction, which would also lead to 
^-compression since the lattice has a positive Poisson ra¬ 
tio set by P). In other domain wall geometries, or other 
orientations of the domain wall relative to the lattice, 
the localized states of self-stress could have small overlap 
with affine strains, which would make them inconsequen¬ 
tial to the buckling behaviour. Alternatively, states of 
self-stress which overlap significantly with shear deforma¬ 
tions could also be realized which would enable buckling 
to be triggered via shear. 


II. BUCKLING IN TOPOLOGICAL CELLULAR 
METAMATERIALS 

The cellular metamaterial differs from the ideal frame 
in two important ways. First, real structures terminate 
at a boundary, and loads on the boundary are no longer 
equilibrated by states of self-stress, but rather by tension 
configurations that are in equilibrium with forces on the 
boundary nodes m- However, these tension states are 
closely related to system-traversing states of self-stress 
in the periodic case, with the forces on the boundary 



FIG. 3. Stretching, shear and bending contributions 
to the linear in-plane response of the cellular meta- 
material. Response of a planar cellular structure, related to 
the lattice in Fig. 2 but with free edges, subject to a vertical 
compressive force F (solid arrows) at each point highlighted 
along the top and bottom edges. The structure is modeled as 
a network of flexible beams connected by rigid joints at the 
nodes, and with each beam providing torsional stiffness in ad¬ 
dition to axial stiffness. The beams are coloured according to: 
a, axial compression; b shear load; c, bending moment. 


nodes in the finite system playing the role of the tensions 
exerted by the boundary-crossing beams in the periodic 
system. Therefore, the states of self-stress also provide 
information about the load-bearing regions in the finite 
structure away from the boundary. 

In addition to having boundaries, the cellular block 
probed in Fig. 1 also departs from the limit of an ideal 
frame, as it is made of flexible beams rigidly connected at 
the nodes and can support external loads through shear 
and bending of the beams in addition to axial stretch¬ 
ing or compression. Nevertheless, the states of self-stress 
and tension states of the corresponding frame (with the 
same beam geometry) determine the relative importance 
of bending to stretching in the load-bearing ability of 
the cellular structure [21]. To verify that the localized 
states of self-stress in the underlying frame influence the 
response of the finite cellular structure, we numerically 
calculate the in-plane response of each layer treated as an 
independent 2D cellular structure with loading confined 
to the 2D plane. Each beam provides not just axial ten- 
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sion/compression resistance but also resistance to shear 
and bending. A beam of length L, cross-sectional area A 
and area moment of inertia I resists (i) axial extensions e 
with a tension t = {EA/L)e^ (ii) transverse deformations 
Cg with a restoring shear force s = V1[E1 and 

(iii) angular deflections of the end nodes 0 with restor¬ 
ing moment m = {EIIL)0 (Fig. 3). Since A ^ w‘^ and 
I ^ for a beam of width re, the relative contribu¬ 
tion of the bending, shear and torsional components to 
the total stiffness is set by the aspect ratio wjL oi each 
beam; the frame limit with no bending or shear stiffness 
is recovered when wjL 0. The physical sample has 
beams with aspect ratios in the range 0.1 < wjL < 0.2, 
indicating a small but appreciable contribution of shear 
and torsion to the response. 

The 2D linear response of such a structure is calculated 
by augmenting the equilibrium matrix A to include an 
additional degree of freedom (a rigid rotation angle) at 
each node, and two additional restoring forces (shear and 
torsion) for each beam, see Materials and Methods. The 
resulting equilibrium matrix is of size Srin x 3nb, which 
implies that cellular solids based on isostatic frames with 
rib ~ drin are severely over constrained for d = 2,3, and 
structurally stable even with free boundaries. In particu¬ 
lar, they can support any loads exerted on the boundary 
nodes as long as the net forces and torques about the 
center of mass are zero. Once the equilibrium matrix 
is constructed, its singular value decomposition can be 
used to completely determine all stresses and torques ex¬ 
perienced by the beams in response to forces specified 
on the boundary nodes m Fig. 3 shows the linear re¬ 
sponse of each plane of the cellular pattern correspond¬ 
ing to the domain wall geometry of Fig. 2 under uniform 
compressive force applied to the nodes on the top and 
bottom edges. The beams participating in the states of 
self-stress at the left domain wall are singled out by their 
high compression, whereas the rest of the structure pri¬ 
marily supports the boundary load through shear rather 
than compression or bending. Remarkably, the unique 
compression-dominated response of the left domain wall 
(originating from the topological invariant in the ideal¬ 
ized isostatic frame with freely-hinged joints) survives in 
cellular structures away from the isostatic condition. 

We expect a similar localized compression-dominated 
response in each layer of the stacked structure (Fig. Ib- 
c). The enhanced compressions along the left domain 
wall trigger buckling when the compression exceeds the 
Euler beam buckling threshold, U < —cEljLl^ where c 
is a positive numerical factor determined by the clamping 
conditions as well as cooperative buckling effects. Buck¬ 
ling is signified by a loss of ability to bear axial loads, 
as the beam releases its compression by bending out of 
plane. Upon compressing the 3D sample between two 
plates as shown in Fig. 4a, we see a significant out-of¬ 
plane deflection for beams along the left domain wall 
(Fig. 4b-c), consistent with buckling of the maximally 
stressed beams in Fig. 3a. The deflection in different 
layers is coordinated by the vertical beams connecting 


equivalent points, so that beams within the same column 
buckle either upwards (column 1) or downwards (column 
2) to produce a distinctive visual signature when viewed 
along the compression axis. Beams connecting different 
planes in the stacked pattern create additional states of 
self-stress that traverse the sample vertically, but these 
do not single out any region of the material, and do not 
couple to the specific in-plane loading of Fig. 4a. 

We emphasize that having as many states of self-stress 
as there are unit cells along the domain wall is crucial for 
the buckling to occur throughout the domain wall. When 
a beam buckles, its contribution to the constraints of dis¬ 
tances between nodes essentially disappears, reducing the 
number of load-bearing configurations by one. If there 
were only a single localized state of self-stress in the sys¬ 
tem, the buckling of a single beam would eliminate this 
state, and the compressions on the other beams would 
be relaxed, preventing further buckling events. However, 
the presence of multiple states of self-stress, guaranteed 
in this case by the topological origin of the states, al¬ 
lows many buckling events along the domain wall. In the 
SI Text, we use an adaptive simulation that sequentially 
removes highly compressed beams to show that each re¬ 
peating unit along the y direction experiences a unique 
buckling event even when the loss of constraints due to 
other buckling events is taken into account (Supplemen¬ 
tary Figure SI). 


III. ROBUSTNESS OF THE BUCKLING 
REGION 

Finally, we show that the robustness of the topological 
states of self-stress predicted within linear elastic theory 
carries over to the buckling response in the non-linear 
regime. There is a wide range of distortions of the de¬ 
formed kagome unit cell which do not close the acoustic 
gap, leaving the topological invariants and hence the 
polarization Rt unchanged HU- These distortions may 
bring the unit cell arbitrarily close to the regular kagome 
lattice with equilateral triangles of beams (for which the 
gap is closed and the polarization is undefined). The 
unit cell of the two dimensional lattice shown in Fig. 5a 
is minimally distorted away from the regular kagome lat¬ 
tice, but this barely noticeable distortion (inset of Fig. 
5a) is sufficient to induce the same topological polariza¬ 
tion Rt in the underlying frame as before. As a result, 
the domain wall on the left (characterized by a net out- 
flux of polarization) localizes states of self-stress (shown 
in Supplementary Figure S2) even though the unit cells 
are nearly identical on either side. 

For ease of visualization, we tested this design in a 
2D prototype cellular metamaterial, obtained by laser¬ 
cutting voids in a 1.5 cm thick slab of polyethylene foam 
(see Materials and Methods), leaving behind beams of 
width 1-2 mm and lengths 10-12 mm (Fig. 5a). The as¬ 
pect ratio of the beams is comparable to that of the 3D 
sample. Since the slab thickness is much larger than the 
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FIG. 4. Buckling in the 3D topological cellular metamaterial. a, Top view of the 3D sample whose construction was 
outlined in Fig. 1, showing the compression applied by confining the sample between transparent plates in contact with the 
front and back surfaces. The buckling zone is highlighted in red as in Fig. 1. Two vertical columns within this zone are labelled 
in yellow. Magenta arrows show the compression direction, b— c, View along the compression axis at compressions of 0 and 
20% respectively. The beams in the region with states of self-stress have buckled in the vertical direction, whereas other beams 
have largely deformed within their stacking planes. Scale bar, 25 mm. 


beam width, deformations are essentially planar and uni¬ 
form through the sample thickness, and can be captured 
by an overhead image of the selectively illuminated top 
surface. The sample was confined between rigid acrylic 
plates in free contact with the top and bottom edges, and 
subjected to a uniaxial in-plane compression along the 
vertical direction by reducing the distance between the 
plates. Using image analysis (see Materials and Meth¬ 
ods) we computed the beams’ tortuosity, defined as the 
ratio of the contour length of each beam to its end-to-end 
distance. Buckled beams have a tortuosity significantly 
above 1. Under a vertical compression of 4%, only beams 
along the left domain wall show a tortuosity above 1.05, 
consistent with a localized buckling response (Fig. 5b). 
However, under higher strains, the response of the lattice 
qualitatively changes. When beams have buckled along 
the entire domain wall, its response to further loading is 
no longer compression-dominated (Supplementary Fig¬ 
ure SI). Future buckling events, which are triggered by 
coupling between torsional and compressional forces on 
beams due to the stiff hinges, no longer single out the 
left domain wall and happen uniformly throughout the 
sample (Fig. 5c and Supplementary Movie 2). 

We have demonstrated that piling up localized states 
of self-stress in a small portion of an otherwise bending- 
dominated cellular metamaterial can induce a local 
propensity for buckling. Whereas this principle is of gen¬ 
eral applicability, our buckling regions exploit topolog¬ 
ical states of self stress m which provide two advan¬ 
tages. First, they are indistinguishable from the rest of 
the structure in terms of node connectivity and material 
parameters, allowing mechanical response to be locally 
modified without changing the thermal, electromagnetic 
or optical properties. This feature could be useful for op¬ 
tomechanical [30] or thermomechanical m metamaterial 
design. Second, the buckling regions are robust against 



FIG. 5. Buckling is robust under polarization¬ 
preserving changes of the unit cell, a, A 2D foam cellular 
prototype, whose unit cell maintains the topological polariza¬ 
tion Rt even though its distortion away from the regular 
kagome lattice is small (a zoom of the constituent triangles is 
shown within the yellow circle). The domain wall geometry is 
identical to that of the 3D sample, with the left domain wall 
localizing states of self-stress. Scale bar, 2 cm. b, Response of 
the structure under a vertical compression of 4% with free left 
and right edges. The beams are coloured by the tortuosity, 
the ratio of the initial length of the beam to the end-to-end 
distance of the deformed segment (colour bar), c, Response 
of the structure under 7% compression, with beams coloured 
by tortuosity using the same colour scale as in b. 
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structural perturbations, as long as the acoustic bulk gap 
of the underlying frame is maintained. This gap is a prop¬ 
erty of the unit cell geometry. Large deformations that 
close the gap could be induced through external actua¬ 
tion [32] or confinement, potentially allowing a tunable 
response from localized to extended buckling in the same 
sample, reminiscent of electrically tunable band gaps in 
topological insulators ESI. 


IV. MATERIALS AND METHODS 

A. Deformed kagome lattice unit cell 


The deformed kagome lattices are obtained by deco¬ 
rating a regular hexagonal lattice, built from the primi¬ 
tive lattice vectors {ai = ax, 3.2 = —(a/ 2 )x + {^/3a/2)y} 
where a is the lattice constant, with a three-point unit cell 
that results in triangles with equal sizes. The unit cells 
are described by a parametrization introduced in Ref. mi 
which uses three numbers (xi, X 2 , X 3 ) to quantify the dis¬ 
tortion of lines of bonds away from a regular kagome 
lattice {xi = 0). The unit cell in Fig. la, which forms 
the basis for the 3D cellular structure, is reproduced by 
(xi,X 2 ,X 3 ) = (—0.1,0.06,0.06) for which the topological 
polarization is Rt = —ai [T3|. The unit cell for the de¬ 
sign in Fig. 5a parametrized by (—0.025,0.025,0.025) has 
the the same polarization. In all cases we describe the 
unit cell and polarization in the outer region; the inner 
region between the domain walls rotates this unit cell by 
TT which flips the polarization direction. 


B. Construction of the equilibrium/compatibility 
matrix 


Analysis of the linear response of a frame or a cellular 
material begins with the construction of the equilibrium 
matrix A. Since it is more natural to relate node dis¬ 
placements to beam length changes, we describe how to 
build the compatibility matrix C = A^, relating point 
displacements u to extensions e via Cu = e. We con¬ 
struct the compatibility matrix from the contributions of 
individual beams. Consider a single beam aligned to the 
X axis connecting hinge 1 at (0,0) to hinge 2 at (T,0). 
There are six degrees of freedom potentially constrained 
by the beam: the positions (x^,^^) at each node i and 
the torsion angle Oi. Within Euler-Bernoulli beam the¬ 
ory, these values are sufficient to determine the shape of 
the beam along its length as well as the forces and torques 
at each end needed to maintain equilibrium [34|. Three 
independent combinations of forces and torques can be 
identified which are proportional to generalized strains 
experienced by the beam: pure tension t oc X 2 — xi, pure 
shear s oc y 2 —yi—L{ 0 i-\- 02 )/ 2 ^ and pure bending torque 
m (X O 2 —O 1 , 8iS illustrated schematically in Fig. 3. In ma¬ 


trix form, this gives 


/-I 0 0 1 0 0 \ 

0 -1 -L/2 0 1 -L/2 

\ 0 0 -1 0 0 1 / 


/Xl\ 

yi 

0i 

X2 


= Cu. 


2/2 

\ 02 / 


(6) 


The forces and torque (generalized stresses) are obtained 
from the generalized strains e through the stiffness ma¬ 
trix K which depends on the Young’s modulus E, the 
cross-sectional area A, the area moment of inertia I and 
the beam length L: 


/t\ lEA/L 0 0 \ 

5 = 0 \2E1/L^ 0 . e = Ke. (7) 

\m) \ 0 0 EljL) 


The compatibility matrix for a beam with arbitrary ori¬ 
entation is obtained by projecting the displacement vec¬ 
tors at the end of each beam onto the axial and transverse 
directions using the appropriate rotation matrix which 
depends on the angle made by the beam with the x axis. 
Each beam in an assembly provides three rows to the 
compatibility matrix, with additional columns set to zero 
for the degrees of freedom unassociated with that beam. 
In the simpler frame limit, each beam only resists axial 
extensions , and contributes one row (the first row of the 
C matrix in Eq. to the compatibility matrix. 

Once the equilibrium matrix is constructed, the ap¬ 
proximate states of self-stress of the periodic frame 
(Eig. 2) as well as the linear response of the cellular ma¬ 
terial under loads (Eig. 3) are obtained from its singular 
value decomposition following the methods of Ref. m 
More details of the computation are provided in the SI 
Text. 


C. Construction and characterization of 2D and 
3D prototypes 

The 3D structures were printed by Materialise N.V. 
(Leuven, BE) through laser sintering of their proprietary 
thermoplastic polyurethane TPU 92A-1 (tensile strength 
27 MPa, density 1.2 g/cm^. Young’s modulus 27 MPa). 

The 2D structures were cut using a VersaLaser 3.5 
laser cutter (Laser & Sign Technology, New South Wales, 
AU) from 1.5 cm thick sheets of closed-cell cross-linked 
polyethylene foam EKI-I306 (EKI B.V., Nijmegen, NL; 
tensile strength 176 kPa, density 0.03 g/cm^. Young’s 
modulus 1.7 MPa). A characteristic load-compression 
curve of a 2 D sample is shown in Supplementary Eigure 
S3. 


D. Image analysis of 2D experiment 

Images of the 2D cellular prototype (Eig. 5) were ob¬ 
tained using a Nikon CoolPix P340 camera and stored 




as 3000x4000 px 24-bit JPEG images. To quantitatively 
identify the buckled beams in the 2D prototype under 
confinement, we extracted the tortuosity r of each beam, 
defined as the ratio of the length of the beam to the dis¬ 
tance between its endpoints. Tortuosity was estimated 
from the sample images through a series of morphologi¬ 
cal operations, as detailed in SI Text and Supplementary 
Figure S4. Straight beams have r = 1 whereas beams 
that buckle under axial compression have r > 1. Un¬ 
like other measures of curvature, tortuosity distinguishes 
buckled beams from sheared beams (schematic. Fig. 3b) 


which have r > 1. 
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Supplementary Information 

I. SVD ANALYSIS OF EQUILIBRIUM MATRIX 

The linear response of the frame as well as the cellular 
metamaterial are obtained using the singular value de¬ 
composition (SVD) of the equilibrium matrix, as detailed 
in Ref. [15] and summarized below. The SVD analysis si¬ 
multaneously handles the spaces of node forces as well as 
beam stresses, and properly takes into account states of 
self-stress in the structure. 

The SVD of the equilibrium matrix A with rij- rows and 
Tie columns expresses it as the product of three matrices: 

A = UVW^, (SI) 

where U = [ui,...,Unr] is a square matrix composed of 
nr independent column vectors u^, W = [wi, ..., is 
a square matrix of ric independent column vectors w^, 
and V is a matrix with r non-negative values vu on the 
leading diagonal and all other values zero, r being the 
rank of A. 

The first r columns of U provide the finite-energy 
displacement modes of the structure, whereas the remain¬ 
ing Ur — r column vectors are the zero modes. Similarly, 
the first r and remaining ric — r column vectors of W 
provide the finite-energy bond tension configurations and 
the states of self-stress respectively. 

For the frame in main text Fig. 2 under periodic bound¬ 
ary conditions, the count of topological states of self¬ 
stress predicts eight localized states of self-stress if the 
left domain wall were isolated, in addition to the two 
extended states of self-stress expected under periodic 
boundary conditions. However, since the separation be¬ 
tween the left and right domain walls is finite, the SVD 
produces only two actual states of self-stress, and six ap¬ 
proximate states of self-stress with large va. The corre¬ 
sponding tension configurations are the last 

eight column vectors of W. For of these (w^^_ 7 , w^^_6, 
and w^c) have a significant overlap with Sxx or 
Syy. These are displayed in Fig. 2a-d respectively. 

The response of the cellular network in Fig. 3 of the 
main text is calculated for the finite block with free edges. 
Unlike the frame, this is highly overconstrained and has 
many generalized states of self-stress (now corresponding 
to mixtures of tensions, shears and torques that maintain 
equilibrium). An external load 1 , corresponding to a force 
specification on each point, is supportable by the struc¬ 
ture provided its overlap with the space of zero-energy 
motions is zero; i.e. 

[Ur+i, = 0 . (S2) 

Since the cellular block only has three zero motions cor¬ 
responding to the two translations and one rotation, any 
set of forces with no net force or torque on the structure, 
such as the force configuration shown in Fig. 3a, satis¬ 
fies this condition. The generalized stresses cr (three per 
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beam) are then given by 


r , 

E Ui • 1 

11 -■ 


Wi + Wn,-rX, 


(S3) 


i=l 


where the second term is an arbitrary combination of the 
states of self-stress given by the last ric—r columns of W. 
The weights x are determined by requiring zero overlap 
of the generalized strains with the incompatible strains 
(states of self-stress), which leads to the equation 

WT = -W^ _^K-1 (S4) 

i=i 

which can be solved for x to get the complete generalized 
stresses (three per beam, plotted separately in main text 
Figs. 3a-c). 


II. SEQUENTIAL REMOVAL OF BEAMS 


Fig. 3 (main text) showed that the axial compressions 
in the cellular material under uniaxial loading were con¬ 
centrated along the left domain wall, matching the ex¬ 
pectation from the analysis of the frame with similar 
beam geometry. However, every time a beam buckles, 
its axial load-bearing ability is lost. This fundamentally 
changes the load-bearing states of the underlying frame, 
and there is no guarantee that the beams with the high¬ 
est compression continue to be along the left domain wall. 
However, the number of topological states of self-stress 
localized to the left domain wall in the underlying peri¬ 
odic frame grows linearly with the length of the domain 
wall m, suggesting that the load-bearing ability of the 
left domain wall may not vanish completely with just a 
single buckling event. Here, we numerically show using 
a sequential analysis of the finite cellular frame that the 
left domain wall can support as many buckling events as 
there are repeating units along the y direction. 

To recap. Fig. 3a (main text) shows the compressions 
in a finite cellular block, obtained by tiling the same pat¬ 
tern three times along the y direction, under forces ap¬ 
plied to undercoordinated points at the top and bottom 
edge. Since the threshold for buckling under compression 
scales as 1/T? for beam i of length we expect that the 
beam with the largest value of —tiL‘^ will buckle first as 
the compressive force F is increased. We examine the in¬ 
fluence of the buckling event on the compression response 
by removing this beam from the cellular structure and 
recalculating the axial tensions t without changing the 
forces on the boundary points. This process can be re¬ 
peated, sequentially removing the beam with the largest 
value of —tiL‘1 after the linear response of the cellular 
structure has been recalculated. The result is show in 


Fig. SI, in which beams are coloured by —i.e. 
the beam with the highest intensity in each iteration is 
removed in the next. For the first three iterations of the 
sequential process, the beam chosen for removal lies along 


the left domain wall (Figs.[ST|i-c), showing that the mul¬ 
tiplicity of load-bearing states localized to the left domain 
wall is sufficient for at least one buckling event to occur 
for each repeating unit of the domain wall. 

Furthermore, after each repeating unit along the y di¬ 
rection has experienced buckling, the beams in the rest 
of the lattice experience much lower compressions rel¬ 
ative to the initial compressions along the left domain 
wall (Fig. [sT]i), signifying that the lattice remains en¬ 
tirely bending/shear-dominated away from the domain 
wall. 

By performing simulations on different system sizes, 
we confirmed that upon increasing the sample size along 
the y direction, the number of buckling events localized 
to the left domain wall increases proportionally. 


III. LOAD-COMPRESSION CURVE OF A 2D 
SAMPLE 


A rudimentary measurement (supplementary fig¬ 


ure S3) of the load-compression curve of a 2D cellular 
prototype has been performed to assess the effect of the 
localized buckling events on the mechanical response. 
The prototype, shown in Fig. |S3)3 , uses a unit cell char¬ 
acterised by (xi,X 2 ,X 3 ) = (—0.085,0.085,0.085) in the 
parametrization of Ref. M- This is geometrically simi¬ 
lar to the unit cell used for the 3D sample in the main 
text and has the same polarization. However, we use 
more unit cells in the 2D sample to reduce edge effects 
on the measurement. The overall geometry is similar to 
that of the finite 2D lattice used in the linear response 
calculations of main text Fig. 3. 

The 2D cellular foam sample (460 mm x 150 mm x 15 
mm; beam widths 1.5-2.5 mm; see Materials and Meth¬ 
ods for more details) was positioned with its narrow edge 
against a supporting acrylic plate and scale (Kern PCB 
10000-1). The sample was confined from above by a sec¬ 
ond acrylic plate. No lateral confinement was provided. 
The distance between top and bottom plates was incre¬ 
mentally decreased from 140 to 90 mm using a laboratory 
jack. At each confinement step the system was allowed 
to relax and the equilibrium force on the sample edge 
was recorded from the scale, resulting in a quasistatic 
measurement of the sample’s load-compression curve. A 
small upward drift of ca. 0.5 kN was recorded in the 
scale’s force response due to long relaxation times. This 
drift (assumed to be linear in the confinement) has been 
subtracted from the recorded data. The complete mea¬ 
surement is shown in figure |S3^ . Buckling of the beams 
along the domain wall leaves a clear signature in the form 
of a sudden softening of the response (region highlighted 
in red). Prior to buckling, the main resistance to con¬ 
finement was provided by the compressed beams along 
the domain wall since the rest of the beams could easily 
bend to accommodate the confinement. Buckling effec¬ 
tively removes the compressional stiffness of the beams 
along the domain wall. Further confinement is only re- 
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FIG. SI. Sequential removal of beams to recreate the effect of buckling in the linear response, a, Axial compression 
of the cellular material shown in Fig. 3a under identical loading, but with beams coloured by the propensity for buckling, 
—tiL^jF. The beam with the highest value of this quantity will be the first to buckle as the force F on the boundary points 
is ramped up. b—d, Result of sequentially removing the beam with the highest propensity for buckling and recalculating the 
compressions in the remaining beams under the same axial loading. 


sisted by additional bending of all beams, which leads to 
a significantly softer response. 

A linear load-compression followed by a plateau is typ¬ 
ical of elastomeric cellular materials m Such materi¬ 
als also tend to have a secondary stiffening regime due 
to densification when adjacent beams contact each other 
along their lengths. In our measurement, however, the 
entire sample buckled in the lateral direction before this 
regime could be reached. 


IV. DETAILS OF 2D IMAGE ANALYSIS 

In the Materials and Methods section of the main text, 
we outlined the quantitative analysis of the 2D exper¬ 
imental data in Fig. 5 (main text). Here, we provide 
more details of the image processing required to measure 
the tortuosity of each beam in a cellular sample under 
compression. The tortuosity is defined as the ratio of the 
length of each beam to its end-to-end distance. 


The analysis starts with a raw image such as in Fig. 5a 
(main text), a top-down view of the sample whose top 
surface is selectively illuminated so that the beams mak¬ 
ing up the cellular material sample are visible against a 
dark background. Fig. [S4^ shows a subset of a different 
sample with similar beam lengths and widths, on which 
we illustrate the image processing steps. 


A. Isolating beam sections 

To quantify the tortuosity of the lattice beams, the 
neutral axis (center line) of each beam must be iso¬ 
lated. The image is first binarized through intensity 
thresholding (figure |S^): pixels with an intensity value 
above/below a certain threshold are set to 1/0. This sep¬ 
arates the structure (1) from the background (0). 

The finite-width beams are then reduced to 1-pixel¬ 
wide center lines by finding the morphological skeleton of 
the binarized image. To do this, we perform a repeated 
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FIG. S2. Localized states of self-stress for frame corresponding to design of Fig. 5a (main text). The numerically- 
obtained states of self-stress for a frame with the same domain wall geometry as in Fig. 2, but with the unit cell parametrized 
by {xi,X 2 ,X 3 ) = (—0.025,0.025,0.025) which has a smaller distortion away from the regular kagome lattice and is used in 
Fig. 5a. The left domain wall still has a net outflux of the topological polarization and localizes states of self-stress, which have 
a significant overlap with the affine bond extensions Cx and Cy associated with uniform strains along x and y respectively. 



Total strain 
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FIG. S3. Measurement of the total confining force as a function of imposed strain for a 2D cellular prototype 
sample (Fig. 5). a, The measured data (black circles) and an estimate of the measurement error (dark grey contour) are 
shown. The red region indicates the strain regime in which the beams within the topologically rigid domain wall buckle. In 
the light grey region, the sample response is dominated by out-of-plane buckling of the entire sample rather than by in-plane 
deformations of the cellular structure. The error was estimated at 0.2 kN based on uncertainties in the compression and force 
due to the rudimentary measurement apparatus, b, 2D cellular prototype for which the load-compression curve was measured. 
Its unit cell in the outer region is parametrized by (xi, X 2 , X 3 ) = (—0.085,0.085,0.085). 


sequence of binary morphological operations: skeletoniz¬ 
ing and spur removal. The former consists of removing 
edge ^ 1 ’ pixels without breaking up contiguously con¬ 
nected regions. The latter is a pruning method that re¬ 
moves protrusions under a certain threshold length; it 
reduces spurious features due to uneven lighting and res¬ 


olution limitations. The resulting skeleton is shown in 
figure [S4|:. 

Next, we break up the skeleton into the neutral axes of 
individual beams. Note that the topology of the frame, 
with four beams coming together at each joint, is not 
reproduced in the skeleton, due to the finite size of the 
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beam joints. Instead, the skeleton has short segments 
connecting junctions of pairs of beams. To remove these 
segments, we identify the branch points of the skeleton 
(pixels with more than two neighbors) and clear a circular 
area around them. The size of the segments and of the 
cleared disk is of the same order as the original beam 
width. The resulting image, consisting of the beams’ 
approximate neutral axes, is shown in figure [S4| 

Once the beam axes have been separated, two steps are 
implemented to improve the signal to noise ratio. First of 
all, beams near the edges of the sample are ignored. They 
are not always separated correctly from their neighboring 
beams, since they have lower connectivity than beams in 
the bulk of the lattice; this can lead to overestimation 
of the tortuosity. Secondly, since our lattices are fairly 
regular, the beams have a restricted size. Beam elements 
above and below these length thresholds are artifacts, 
and we ignore them. Figure [S^ illustrates the processed 
image data after these final steps. 


B. Computing the tortuosity of each neutral axis 

After the beam axes have been isolated, connected 
components (sets of contiguous ^ 1 ’ pixels) are found and 
labeled. Each connected component corresponds to an 
individual beam axis. Numerical data can be extracted 
from it by analyzing its pixels’ coordinates, rj = (xp^yp). 
These coordinates are fitted with polynomials (x (/), ^(/)) 
parametrized by the normalized path length along the 
beam, I G [0,1]: 



Cubic polynomials in I successfully capture the shape of 
the neutral axis of each beam, while smoothing out pixel 
noise which leads to small-wavelength features in the bi¬ 
nary skeleton. Fig. [S4] F displays the smoothed contours 
superimposed on the original image, showing that the 
shapes of the individual beam sections are faithfully re¬ 
produced by the parametric form of Eq. [85] 

The tortuosity r of the beam axes is the arc length 
A along the line divided by the distance D between the 
endpoints. Both quantities can be extracted from the 
fitted model of the beam axes, f{l) = {x{l),y{l)): 



D=|r(l)-f(0)| (S7) 


Once the tortuosities of the beam segments have been 
computed, pixels in each segment of the binary image 
(EigurejS^) are coloured according to the corresponding 
tortuosity to obtain images such as Eigs. 5b-c of the main 
text. 
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(a) Example of experimental data 
image (cropped). 



(b) The image is binarized through 
intensity thresholding. 



ZL Zi_iZ JZ 


(d) Node segments are removed. 


(e) Image edges and noise features 
are cleared. 



(c) The image is skeletonized 
through morphological operations. 



(f) The fitted curves are shown, su¬ 
perimposed on the original image 
and colored according to their tor¬ 
tuosity. 


FIG. S4. Demonstration of steps in the image analysis process. 
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buckles out of the layer producing a distinctive visual 
signature. 



Supplementary Movie 1: View of the 3D layered cel¬ 
lular metamaterial prototype as it is compressed into the 
plane of the figure*. A subset of beams singled out by 
topological states of self-stress in the underlying frame 


Supplementary Movie 2: View of the 2D cellular 
metamaterial prototype as it is compressed vertically. 
Along with each raw image, a reconstructed image with 
beams coloured by their tortuosity (as in main text 
Fig. 5) is also shown. 













